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Variational study of Fermi-surface deformations in Hubbard mod- 
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Abstract -We study the correlation-induced deformation of Fermi surfaces by means of a new 
diagrammatic method which allows for the analytical evaluation of Gutzwiller wave functions in 
finite dimensions. In agreement with renormalization-group results we find Pomeranchuk insta- 
bilities in two-dimensional Hubbard models for sufficiently large Coulomb interactions. 



Introduction. — The shape of the Fermi surface (FS) 
is one of the most important properties that determine the 
low-energy physics of electron liquids. The single-particle 
energy levels of non-interacting electrons depend on the 
crystal momentum k from the Brillouin zone through the 
(single-band) dispersion relation e(k). For N electrons at 
zero temperature, all single-particle states which lie below 
the Fermi energy E-p are occupied. The FS separates occu- 
pied and unoccupied single-particle levels, i.e., it consists 
of all k-points which obey the equation Ep — e(kp) = 0. In 
the presence of finite Coulomb interactions, the calculation 
of the FS requires the real part of the proper self-energy 
S(k, u) so that the FS is obtained from 



E F - e(k F ) - <Ke£(k F , E F ) = 



(1) 



Within perturbation theory for the Coulomb interac- 
tion PQ, the proper self-energy is defined as the sum over 
all irreducible diagrams for the single-particle Green func- 
tion. 

The calculation of S(k, oS) is a notoriously difficult task 
in correlated-electron theory, even for a single-band Hub- 
bard model 

H = H + h , H Q = ^2 *u4,<r g j,<7 ' ^ = "U™U 

i hi," 

(2) 

in two dimensions. Here, i = (11,12) denotes one of the 
L sites on a square lattice, and a =t, i- The dispersion 
relation is given by 



For small U, the self-energy S(k, u>) may be calculated 
from straightforward perturbation theory (2H1], or using 
renormalisation-group methods [HE]. When U is of the 
order of the bare bandwidth W, or larger, only purely 
numerical methods such as Quantum-Monte Carlo [7J[S] 
and Exact Diagonalisation [9] are available which still suf- 
fer from serious finite-size limitations [9lll0j. In view of 
these significant problems on the theoretical side, our un- 
derstanding of experiments on two-dimensional Fermi sur- 
faces, e.g., those of doped cuprates, is far from satisfactory. 
Moreover, only reliable many-particle approaches permit a 
meaningful comparison of measured and calculated Fermi 
surfaces that may reveal the correct form of Ho in ([2]) ; for 
a recent overview, see, e.g., Ref. [9]. 

In this work, we introduce a new analytical scheme 
to evaluate expectation values for Gutzwiller wave func- 
tions in finite spatial dimensions in a controlled way. By 
construction, the Gutzwiller approach provides the Fermi 
surface of quasi-particles in Landau's Fermi-liquid theory. 
Therefore, the Gutzwiller wave function is an appropri- 
ate tool for the calculation of correlation-induced FS de- 
formations at moderate interaction strengths, U ~ W. 
Unlike numerical schemes for the evaluation of Gutzwil- 
ler wave functions pD - [15] our approach does not suffer 
from finite-size limitations. It therefore provides us with a 
momentum-space resolution that is needed for the study 
of Fermi-surface deformations. 

Evaluation of Gutzwiller wave functions. The 

variational wave functions introduced by Gutzwiller |16j 
for the single-band Hubbard model ([2]) have the form 



(3) 



) =P G |*0) =n.A;G|*0) , 



(4) 
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where \^o) is a (normalised) single-particle product state 
and the local 'Gutzwiller correlator' is defined by 



A G = i - (i - <?M 



(5) 



Here, g > is a variational parameter which allows for 
the optimisation of the average number of doubly-occupied 
lattice sites. In this work, we will use the more convenient 
definition 



Pi 



r 



A r |r>ii<r| 



(6) 



for the local correlator. It contains the variational param- 
eters Ar for the four local states 



|r> i e{|0) i ,|t)i,li)i,lti)i} 



(7) 



for the empty, singly, or doubly occupied site i. To keep 
notations simple, we assume that the wave functions |^g) 
and |\&o) are translationally invariant and paramagnetic. 
The corresponding derivation for more general wave func- 
tions is straightforward. 

In order to determine the expectation value of the 
Hamiltonian (|5J) with respect to the wave function (j4]) we 
need to evaluate (i ^ j) 



(*g|*g> 
(*g|*|*g> 



<*giCaj*g> 



pAPi n p i 



n p < 



l(#iJ) 



(8) 
(9) 
(10) 



where (. . .)o denotes expectation values with respect to 



|*o) and c{ 



As we will show below, our 



diagrammatic expansion significantly simplifies if we set 



where 



af F 



Pi 



-HF ~HF 



1 + xd 



HF 



n l,t n U ' 



U l.cr 



(11) 



(12) 



and no = (ni,<r)o = N/(2L). Equation determines 
three of the four parameters Ap as well as the coefficient x. 
In this way, we are left with only one variational param- 
eter. For instance, we may express the parameters A0, 
Ai = A CT for empty and singly-occupied sites by A^, the 
parameter for doubly-occupied sites. Alternatively, due 
to the relations 

x = [Xl ~ !]/(! - «o) 2 & \ 2 d = 1 + x(l - n Q ) 2 , (13) 

we may also consider x as our variational parameter. The 
expansion (fTTj) was first introduced in Ref. [17] where it 
has been used as a convenient tool for the evaluation of 
expectation values in infinite dimensions D = 00 and of 
1/D corrections. For another series expansion for the 
Gutzwiller wave function around the limit D = 00, see 
Ref. [T8] , Note that both local correlators Pi-G and Pi, 



together with the condition (fTTj) . lead to the same vari- 
ational space if the single-particle wave function j^o) is 
treated variational object. 

With Eq. (TTTT) . the norm and the expectation values 
Eqs. (|5j)- (|10[) are given in form of a power series in x, 

00 k 

<*Gl*G> = £^E'< d r..aJo' (14) 

fc=0 li,...l fe 
00 

<*G|^|*G> = ^EFE'<^...I fc >0. ^ 
fe=0 'h,...l k 

00 fa 

(* g \4a\* g ) - Ef E'( ? U^ f ..,iJ '( 16 ) 

k=0 1 ....1. 



where we introduced the notation 



iHF = ?HF . , , iHF 

- "ll "lfc 



jHF _ 



1 



(17) 



The primes in Eqs. (|14[) - (|16[) indicate the summation re- 
strictions l p 7^ l p >, l p ^ i,j for all p,p'- Note that the 
same expressions arise for the original Gutzwiller correla- 
tor Pi : G when we replace x and d^ F by (g 2 — 1) and d\, 
respectively |19|I20|. As we will demonstrate below, our 
expansion with respect to x converges significantly faster 
than the expansion in (g 2 — 1). Therefore, the first few or- 
ders in the x-expansion permit us to evaluate the Gutzwil- 
ler wave function accurately for not too large interaction 
strengths. 

The expectation values in Eqs. (TT4")) - (|16p can be eval- 
uated by means of Wick's theorem p] . By construction, 
we eliminated all diagrams with local 'Hartree bubbles' at 
internal vertices. To achieve the same for the external ver- 
tices we rewrite the corresponding operators in Eqs. (|15l) . 
OH) as 



di 



(1 - xdo)dr + n (n^ + ) + d i? ,(18) 



42 



where we introduced 

do = 1 

q 



n , 

Ai(A d n + A (l 
Ai(Ad — A0) , 



no)) , 



(19) 



(20) 
(21) 
(22) 



and | =1, I =t- When inserted into ([15]) , the last term 
in ([IS]) combines to A^q^oI^g) so that it does not have 
to be evaluated diagrammatically. In the resulting di- 
agrammatic expansion of Eqs. p^|) - p^)) . the fcth-order 
terms correspond to diagrams with k 'internal' vertices 
on sites li, . . . , U, one (two) 'external' vertices on site i (i 
and j) and lines 



P, 



l.i' 



At 



r )o - <5i,i' n o 



(23) 



connecting these vertices. 

As the final analytical step of our derivation, we ap- 
ply the linked-cluster theorem pQ. The norm ([i"4]) cancels 
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K 3 = '■ ^<S^> J + 
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Fig. 1: Diagrams with up to two internal vertices for /' 4 ' and 
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the disconnected diagrams in the two denominators ([15]) 
and . Note that a straightforward application of this 
theorem is hampered by the summation restrictions in 
these equation. However, applying Wick's theorem to the 
expectation values in (jT4 j) — ffl6 ]l is equivalent to the evalu- 
ation of determinants such as \P\\i | with 1, 1' G (li, . . . , 1^). 
Since these determinants vanish if any of the lattice sites 
li, . . . , lft are the same, we can lift the summation restric- 
tions in (TT4| -- (fT6f without creating additional terms. 

The remaining task is to evaluate the six diagrammatic 
sums 



„k 



s = E ij s W 



fe=0 



(24) 



with 



S E {/( 2 ),/( 4 ),T (1) ' (1) ,T (1) ' (3) ,T (3) ' (1) ,T (3) ' (3) } (25) 



and 



27 



(1)[(3)],(1)[(3)] 



•J 



(*) - E WWJs 



i i. 



Here, (. . . }q indicates that only connected diagrams are to 
be kept. As examples, we show the leading order diagrams 
of and r. ( P' (3) in Fig. [TJ 

The variational ground-state energy functional is given 

by 

(H) G =E Q (\y ),x)=L(E kin + Ud) (27) 



in the form 



(H) G = 2£*y(ff a rS 



2^(1), (1) 



+LUX 2 d ({l - xd )I {4} + 2n I {2) + d ) , (28) 

where |\&o) enters the energy expression solely through the 
lines P?- and through no. In our case of a translationally 
invariant wave function, wc have 



(n ij(T )c 



n 



(29) 



Fig. 2: Difference between the exact double occupancy d c and 
its n-th order Taylor expansions d n (with n = 3,5,7,9, 11 in 
one dimension (in descending order) with expansion parame- 
ters x (solid line) and (g 2 — 1) (dashed line). 



The l.h.s. of (|2"9"|) is given diagrammatically as 
(rii,a>G - X 2 d (d +I ( - 4 \l-xd ) + 2n I^) 



(30) 



+X 2 (m° + 1^ (i _ 2no) - / (4) (1 + x)m\) , 
where m? = ?i (l - n ). From ([29]), ([30]), (JTSJ), and 

Ai = 1 - ra (1 - n )x , (31) 
we then find the following relation between 1^ and I^ 2 \ 



j(2) = _j(4) 



.t(1 — 2no) 
1 + a;no(l - no) 



(32) 



Therefore, only 1^ in ([2"5]) needs to be calculated. 



The one-dimensional Hubbard model. We have 
tested the quality of our x-expansion against the exact 
results for the Gutzwiller wave function in one dimen- 
sion [HIEO]. From the analytical results one can, e.g., 
determine the Taylor expansions of the average double 
occupancy with respect to x or (g 2 — 1) [21) . It turns 
out that the x-expansion converges much faster than the 
expansion in (g 2 — 1). This can be seen, for example, in 
Fig. [5] where the differences between the exact double oc- 
cupancy (d e ) and its n-th order Taylor expansions (d n ) at 
half band- filling (no = 1/2) are shown for both parameters 
as a function of d e . The figure reveals that the cc-expansion 
to third order is by an order of magnitude closer to the 
exact result than the llth-order expansion in (g 2 — 1). 

In order to assess the absolute quality of our diagram- 
matic x-expansion, we prefer to define the 'order' of the 
expansion by the number of internal vertices which we 
retain in the diagrams. Consequently, the corresponding 
expressions for the kinetic energy E^ m and the double oc- 
cupancy dk to /cth order are not identical to the fcth order 
Taylor expansions because the coefficients q, a, Xd also 
depend on x. In Fig. [3] we show the deviations of the 
results to leading orders (k < 4) for the kinetic energy 



p-3 



J. Biincmann, T. Schickling, and F. Gebhard 




Fig. 3: Differences between the exact double occupancy d e 
(right axis) and kinetic energy E^ ln (left axis) within the 
Gutzwiller wave function in one dimension at half band- 
filling [191120] and the corresponding fcth order diagrammatic 
results for the double occupancy dk (solid lines, n — 0, 1, 3, in 
descending order) and the kinetic energy E^. ln (dashed lines, 
k = 0, 2, 4, in descending order). 

and the double occupancy from the analytic results in one 
dimension at half band-filling (no = 1/2) as a function 
of x (x = —4 corresponds to the atomic limit d = 0). 
At half band-filling the lowest-order results (k = 0) are 
equivalent to the Gutzwiller approximation |17(I19| be- 
cause a(no = 1/2) = 0. Moreover, the even (odd) orders 
k > vanish for the double occupancy (kinetic energy). 
Fig. [3] shows that the 4th-order results reproduce the exact 
results very well. 

The one-dimensional Hubbard model is the worst case 
for our formalism because the latter is exact for the Gutz- 
willer wave function in the opposite limit of infinite di- 
mensions D = oo where all diagrams vanish. Therefore, 
we expect that the 4th-order results for the Gutzwiller 
wave function in two dimensions, which we discuss in the 
following, are also quite accurate. 

Fermi-surface deformations in two-dimensional 
Hubbard models. — For the two-dimensional Hubbard 
model, we treat l^o) an d its FS as variational quantities. 
The minimisation of ([28")) with respect to |^o) leads to the 
effective single-particle equation 

^o ff |*o) = ^^c jiff |fo)=F ff |*o), (33) 

s dEo(\* ),x) 

The effective dispersion relation 

e cff (k) = i^cxp i ( i -j) k t?? (35) 

defines the quasi-particle FS via 

e eff (k F ) = E F (36) 




Fig. 4: Polar plot of the FS deformations 8kF for U/t = 10 
and various band fillings 2no; Inset: FS deformations 8k f for 
2no = 0.9 and various values of U ft. 

because the correlated momentum distribution 

n k ,a = (4,<A,<7>G (37) 

has step discontinuities exactly at the momenta given by 
Eq. ((311). 

The remaining problem is to solve self-consistently the 
closed set of equations (|2?[). ([25)1. ([23) -([Ml), together with 

|-£ (|*o),x)=0. (38) 

Note that E is just an auxiliary quantity, the ground- 
state energy of the effective single-particle Hamiltonian 
i?Q ff . It must not be confused with the variational ground- 
state energy (|28[). Numerically, we determine |^o) by solv- 
ing (|33[) in momentum space while the diagrammatic sums 
in ([25)) and the derivatives in (|34[) are evaluated in posi- 
tion space (up to 4th order in this work). For example, in 
a paramagnetic state, the first two contributions of I A in 
Fig. Q] are given by 

= *E p 4 +4E p L p LKm + ■■■ (39) 

where the lines (|23[) were assumed to be spin independent, 
Pi j = Pfj. Note that the combinatorial factor 4 in (j3T))) 
results from the different possibilities to label the lines 
with spin indices a . We have determined these factors by 
means of a computer algorithm because their calculation 
becomes quite involved for higher-order diagrams. In this 
context, it was particularly helpful to use the exact results 
in one dimension to eliminate programming errors. 

In principle, the effective hopping parameters i?? = 
t°xY have to be calculated for arbitrary distances X = 
i\ — ji, Y = %2 — j2 of the lattice vectors i = 1%) and 
j = However, since the calculation of the deriva- 

tives in ([M)) is numerically quite expensive, we take into 
account only the dominant parameters t c ^ Y - bi our cal- 
culations we include the seven hopping parameters with 
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Fig. 5: Pomeranchuk FS for the Hubbard model with nearest- 
neighbour hopping ti t o = — t for U/t = 10 (solid) in comparison 
with the FS for U/t = (dashed) for 2n = 1.0, 0.95, 0.9 (from 
right to left). Inset: energy gain AE in Kelvin (t — 1.0 eV) in 
the symmetry-broken phase as a function of the filling 2no for 
U/t = 10. 

X 2 + Y 2 < 10. To be consistent, we also set Pf Y = 
for all X 2 + Y 2 > 10 in the energy functional Eq(\^o), x). 
This simplifies the numerical evaluation of the diagrams 
and it ensures that the self-consistency equations lead 
to a (local) minimum of Eq(\^q), x). Note that the nu- 
merical calculation of all diagrammatic sums in (|28l) takes 
less than a minute on a present-day desktop computer. 
Once evaluated, the variational ground-state energy can 
be calculated immediately for any value of U . This illus- 
trates that the study of more complicated wave functions 
or multi-band models in two (or even three) dimensions is 
a feasible problem within our approach. 

We first consider a Hubbard model with only nearest- 
neighbour hopping tifi = — t on a square lattice. Since 
in this case the FS deformations 8k p = kp(U) — kp(0) 
are rather small we plot them in a polar diagram as a 
function of the angle $ = arctan (k y /k x ). Due to the 
symmetry of the lattice, we only need to consider angles 
4> € (0,7r/4). In Fig. @] we show 8k,F as a function of $ 
for various hole-dopings 2no < 1 near half filling. Due 
to particle-hole symmetry, it is sufficient to study hole 
dopings and the FS is unchanged at half band-filling. At 
finite doping, the overall feature of the FS deformation has 
a maximum for 2no ~ 0.95 and it becomes negligible for 
densities 2no < 0.75. Note that the ^-dependence of SkF 
is a non-trivial function of the doping. For example, the 
two curves with 2no = 0.9 and 2no = 0.975 have almost 
the same deformation at $ = (i.e., in [0,1] direction) 



Fig. 6: Pomeranchuk FS for the Hubbard model with nearest- 
neighbour hopping ti,o — —t and second-neighbour hopping 
ti.i — —t' = 0.25t for U/t = 10 (solid) in comparison with the 
FS for U/t = (dashed) for 2n = 0.75, . . . , 1.0(0.05) (from left 
to right). Inset: FS for U/t = 10 (solid) and U = (dashed) for 
2no = 0.9, 1.0, 1.1 (from left to right) in the lattice-symmetric 
phase. 

but differ significantly for finite <F In contrast, for fixed 
no, only the slope of the curves becomes smaller when we 
decrease U; see the inset of Fig. |4j 

For sufficiently large values of U/t and close to half 
band-filling, these Fermi surfaces are unstable against the 
Pomeranchuk effect [2"2"H2"5] . i.e., we find minima in the 
variational space which break the discrete C4 symmetry 
of the lattice, in agreement with related numerical work 
on t-J models [26H28] . In Fig. [5] we show the FS for 
2n = 0.9,0.95,1.0 for U/t = 10 (Pomeranchuk phases) 
and U/t = 0. Around half band-filling, |2n - 1| < 0.12, 
the states with broken symmetry are stable. As seen from 
the inset of Fig. [5l the maximal energy gain due to the 
symmetry breaking is of a fraction of room temperature, 
AE/k B w 100 K for t = 1.0 eV. 

The asymmetry of the Pomeranchuk FS at finite doping 
is quite remarkable. The two intersection points, denoted 
as k\ and k 2 in Fig. obey k\ = 7r — k 2 only at half band- 
filling. The Pomeranchuk minima are not continuously 
connected to those without broken symmetries. Therefore, 
in our approach, the transitions between such states as a 
function of U jt are discontinuous, in general. 

Second, we consider the Hubbard model with an addi- 
tional second-neighbour hopping t\.\ = —t' = Q.25t. Even 
in the absence of symmetry breaking, the FS deformations 
are much larger than in the Hubbard model with nearest- 
neighbour hopping only, as seen from the inset of Fig. [B] 
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where we show the lattice-symmetric FS for U/t = 10 and 
U/t = near half band-filling. The Pomeranchuk insta- 
bilities occur mainly for hole doping, 0.70 < 2no < 1.03 
for U/t — 10, see the inset of Fig. \5\ In Fig. [51 we show 
the Pomeranchuk FS for U/t — 10 and the corresponding 
FS for U/t = for densities 0.75 < 2n < 1. 

Note that the results presented in Figs. |4]|6] are cer- 
tainly altered by the inclusion of other symmetry-broken 
phases with, e.g., magnetic or superconducting order. 
Such phases will be investigated in future studies. One 
should also keep in mind that our approach only allows us 
to study states which are adiabatically connected to some 
non-interacting reference system ('Fermi liquids'). This 
excludes, in particular, the investigation of Mott- Hubbard 
insulators. At zero temperature, however, such insulators 
are usually found only in theoretical studies which delib- 
erately neglect ordered states, such as antiferromagnetic 
spin waves. 

Conclusions. — In summary, we have introduced 
a novel scheme for the evaluation of Gutzwillcr wave 
functions in finite dimensions. As a first application, 
we described the correlation-induced Fermi surface defor- 
mations in two-dimensional Hubbard models within the 
Gutzwiller approach. To the best of our knowledge, there 
exists no alternative method to calculate quasi-particle 
Fermi surfaces for interaction parameters U « W in the 
Hubbard model. Our approach can be extended in var- 
ious directions: to study magnetic and superconducting 
order in single and multi-band Hubbard models [29], and 
to calculate dynamical response functions [30H32] . 
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